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Crystallization Close to the Glass Transition: Dynamic heterogeneities 
do not precede crystallization^ 

Sven Dorosz,*** and Tanja Schilling 



o 

(N 

0) 
tin 

in 



o 



Received Xth XXXXXXXXXX 20XX, Accepted Xth XXXXXXXXX 20XX 
First published on the web Xth XXXXXXXXXX 200X 
DOI: 10.1039/bOOOOOOx 

\ In this article, we address the question whether a crystallization event can be predicted based on observations of the mobil- 
ity distribution in a supersaturated melt. We have carried out computer simulations of overcompressed suspensions of hard 
monodisperse ellipsoids and observed their crystallization dynamics. The system was compressed very rapidly in order to reach 
the regime of slow, glass-like dynamics. We find that, although particle dynamics become sub-diffusive and the intermediate 
scattering function clearly develops a shoulder, crystallization proceeds via the usual scenario: nucleation and growth for small 
supersaturations, spinodal decomposition for large super saturations. 

In particular, we compared the mobility of the particles in those regions where crystallization set in with the mobility in the rest 
of the system. We did not find any signature in the dynamics that pointed towards the imminent crystallization event. Dynamic 
[heterogeneity did not precede crystallization. 
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.1 Introduction 

A liquid that is cooled or overcompressed beyond its freezing 
point either crystallizes or forms a glass. If the degree of un- 
■dercooling (or overcompression, respectively) is small and the 
quench is fast, one usually observes crystallization via nucle- 
ation and growth. If the system is quenched more strongly and 
[beyond its spinodal, it crystallizes immediately. These crystal- 
lization mechanisms compete with another possible process, 
the glass transition, which occurs in many materials if they 
■are quenched sufficiently fast. And between these extremes, 
"mixed routes to crystallization appear in general. 
" A supersaturated melt close to the glass transition resembles 
|a liquid in structure, but it differs from a liquid in its dynamical 
behaviour. In particular, on approach to the glass transition a 
^elt diplays "dynamic heterogeneity", i.e. spatial fluctuations 
■in its local dynamical behaviour. 

In this article we would like to address the question whether 
dynamic heterogeneity and crystallization are linked. One 
could assume that regions of correlated mobility are most 
likely to crystallize and hence attempt to predict crystalliza- 
tion sites and times from the spatio temporal structure of the 
melt. We have tried to do this and did not detect any dynamic 
signal that preceeded crystallization. 

The interplay between crystallization and the glass transi- 
tion in colloids has been addressed by several groups in the 
past, but no definite conclusion has been drawn. Before pre- 
senting our results we give a brief overview over related work: 
The first experiments on the crystallization process in suspen- 
sions of colloidal hard spheres were carried out by Pusey and 



van Megen in the 80s and 90s . During the crystallization 
process the authors detected two-step relaxation in the dy- 
namic scattering function, i.e. dynamical properties of a glass. 
This observation stimulated a series of theoretical and exper- 
imental studies on the interplay of crystallization and glassy 
dynamics in colloidal systems. 

Zaccarelli et al. and Pusey et al.i*^ performed computer sim- 
ulations of hard sphere suspensions for various degrees of 
polydispersity and volume fractions T]. For r\ > 0.58, they 
showed that for moderate degrees of polydispersity the over- 
compressed melt displays signatures of a glass and neverthe- 
less crystallization occurs spontaneously. Particles are dis- 
placed by less than a sphere diameter before crystallization 
sets in. The authors concluded that ergodicity is not essen- 
tial for crystallization, contrary to what had previously been 
assumed ^ Furthermore, the authors conjectured three dif- 
ferent regimes of crystallization. For small overcompres- 
sions, T] < 0.54, crystallization proceeds via nucleation, for 
0.54 < T] < 0.58, the spinodal determines the crystallization 
process, and for r\ > 0.58 glassy dynamics affects crystalliza- 
tion and a new modified mechanism is postulated. 
These studies were continued for mono-disperse hard spheres 
by Valeriani et al.^. In particular the influence of the initial 
preparation process was explicitly considered. It was con- 
firmed that, at high volume fractions, particles move less than 
a sphere diameter before crystallization sets in. But in con- 
trast to ref . ^ Valeriani et al. propose to distinguish only two 
dynamic regimes: Nucleation at volume fractions T] < 0.56 
and crystallization via spinodal decomposition at higher vol- 
ume fractions. A modified crystallization process for T] > 0.58 
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was not found. 

In a follow-up paper Sanz et al. ^ investigated the conjectured 
modified crystallization process for systems of monodisperse 
hard spheres. In contrast to ref. ^, a positive feedback pro- 
cess was identified, which drives the growth of the crystal for 
7] > 0.58. At the interface between the glass and the crystal 
the free volume is increased and particles possess a higher mo- 
bility. The process is described by the authors "as a sequence 
of stochastic micronucleation events correlated in space by 
emergent dynamic heterogeneities''. 

Based on similar results obtained for a molecular system^, 
Tanaka and collaborators investigated the crystallization pro- 
cess in various glassy systems by means of computer simu- 
lations and experiments: Konishi et al.^ reported on experi- 
ments with o-terphenyl (OTP). An increased mobility at the 
crystal/glass interface was observed below Tg because of vol- 
ume contraction within the crystal, but no feedback process 
was observed. Kawasaki et al.— presented simulation results 
for polydisperse hard sphere systems. They infer that the tran- 
sition from the colloidal glass to the crystal shows dynamic 
heterogeneities which are of thermodynamic rather than of ki- 
netic origin. 

In summary, although there has been a number of studies on 
the crystallization process in colloids close to the glass transi- 
tion over the past years, no clear picture has emerged yet. 

2 Simulation technique and data analysis 

We have picked mono-disperse hard ellipsoids as a model sys- 
tem, because they have been shown to exhibit a glass tran- 
sition— and their equilibrium phase diagram has been thor- 
oughly characterized— -J^. This is an advantage over most of 
the model systems used in the studies listed above, for which 
these conditions are not fullfilled at the same time. Also, it 
will probably be possible to validate our computational results 
experimentally, as the interest in the experimental realization 
of ellipsoidal suspensions has been growing recently— 

The properties of a suspension of hard ellipsoids of revolu- 
tion (spheroids) are determined by the volume fraction t] and 
the ellispoid aspect ratio. We will focus on a moderate aspect 
ratio for which the system possesses a first order phase transi- 
tion from the fluid phase to a rotator phase, i.e. to a phase with 
crystalline order in particle positions but without orientational 
order of the axes. 

We carried out Monte Carlo (MC) simulations of suspen- 
sions of monodisperse hard ellipsoids at constant number of 
particles N and constant external pressure P (the temperature 
T is constant, too, but as the system is athermal T only enters 
the discussion as a trivial factor). Particles were propagated by 
local translation and rotation moves. The maximum displace- 
ment per MC step was set to 0.03 particle diameters, the maxi- 
mum rotation of the particle axis to 1 .8°. For such small steps. 



the "dynamics" of the system is similar to Brownian dynam- 
ics and produces the same behaviour on long time-scales22r22. 
The system consisted of A/^ = 10386 prolate hard ellipsoids 
with an aspect ratio of a/Z? = 1.25, where a is the length of 
the axis of symmetry and h is the length of any axis in the 
perpendicular plane. To simplify notation, we introduce the 
dimensionless pressure = (where kg is Boltzmann's 

constant). The coexistence pressure is = 14.34 and the co- 
existence volume fraction is 7] = 0.515 
We studied suspensions at constant external pressure = 
21 ... 30, = 40 and = 50. The corresponding chemi- 
cal potential differences between the supersaturated melt and 
the stable crystalline phase ranges from |A/i| = 0.57 to 

|A^| = 1.08^^^. In FIG. Owe illustrate the state points 
studied here (diamonds) within simulation data for (a) the 
equation of state for our system as well as (b) the chemical 
potential difference between the overcompressed melt and the 
stable crystal. 
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Fig. 1 (a) Equation of state for prolate ellipsoids, aspect ratio 
a/b = 1.25. (b) Chemical potential difference A/i between the 
overcompressed melt and the stable crystal. The diamonds indicate 
the supersaturations for which we studied the crystallization process. 

During the simulation we monitored the volume fraction r] 
and the local ^6^6 -bond order parameteri^Si^l. Two ellipsoids 
/ and j were counted as neighbours if their distance r^y < 1.5^? 
and as "bonded" within a crystalline region if q(y{i) • ^6*(j) > 
0.7. nij{i) is the number of bonded neighbours of the /th par- 
ticle. If n}y{i) > 8 we labelled particle / as "crystalline". We 
computed the distributions of q^{i) • ^6* (j) and of niy{i) for the 
bulk liquid, the bulk crystal and for crystallites embedded in a 
liquid to verify that this criterion, albeit isotropic, works in a 
solution of anisotropic particles. Due to the moderate aspect 
ratio and the fact that we quench into the rotator phase, the 
q6q6-method sufficed to detect the crystallites. 

We started the simulation runs from a meta- stable liquid 
configuration at a volume fraction t] = 0.567. At this vol- 
ume fraction, the induction time for crystal nucleation is very 
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long compared to the self-diffusion time of an ellipsoid, hence 
we could "equilibrate" the system in the meta-stable overcom- 
pressed melt state At ^ = 0, the external pressure was in- 
stantaneously set to a higher value. The system responded and 
slowly relaxed to a higher density. After ^10^ MC steps, the 
volume fraction reached a final plateau value for all simulation 
runs that we present here. 



mean square displacements differs from a Gaussian distribu- 
tion. It quantifies the presence and strength of dynamic het- 
erogeneities. For high overcompression we observe that the 
maximum of the NPG is growing, i.e. that dynamic hetero- 
geneities become more pronounced. In agreement with the 
work of Pfleiderer ^ we conclude that the dynamics is glass- 
like for this regime of overcompression. 



3 Results 

3.1 Dynamics of the melt 

After the systems had relaxed to a plateau value in density, we 
computed the average mean-square displacement of a particle 
in the melt and the dynamic structure factor (see FIG. EJ. We 
only took particles into account that were in a liquid environ- 
ment, crystalline particles did not contribute to the data shown 
here. 
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Fig. 2 (a) Mean- square displacement (Ar^) as a function of 
simulation time for different values of pressure P*. In addition a 
linear growth law and a sub diffusive power law are plotted for 
comparison, (b) Dynamic structure factor as a function of 
simulation time for different values of P* . ^max corresponds to the 
first maximum of the static structure factor S{q). 

From the data shown in FIG. |2] we infer that the simulation 
runs at lower values of pressure < 30 are "equilibrated" in 
the meta-stable liquid basin: The mean square displacement is 
linear as a function of simulation time and relaxation times are 
short in comparison to the induction time for crystal growth 
(see FIG. IH discussion follows below). In this regime parti- 
cles diffuse several times their diameter before crystallization 
sets in. Runs at > 40 converge to a plateau in volume frac- 
tion but the system is far from equilibrium. The mean-square 
displacement is sub-diffusive and the dynamic structure factor 
displays slow relaxation. 

In FIG. [3] we show the non Gaussian parameter (NGP) as 
a function of time for individual simulation runs of various 
values of P*. The NGP measures how far the distribution of 
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Fig. 3 The time evolution of single trajectories is analyzed with 
respect to the non Gaussian parameter. The maximum is increasing 
for larger P* and the associated typical time is growing. 



3.2 Crystallization process 
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Fig. 4 Time evolution of the total number of crystalline particles 
for different external pressures P*. 

The time evolution of the total number of crystalline parti- 
cles is presented in FIG.|4]for several typical simulation runs. 
In the case of = 27, there is a long induction time after 
which the total number of crystalline particles grows rapidly. 
Here we are clearly dealing with nucleation and growth. 

In the case of = 40 however, particles diffuse only a 
fraction of their diameter before crystallization sets in. We 
observe a decrease in the growth rate as compared to = 27. 
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At = 50, the free volume per particle is too small to allow 
for successful rearrangements on the time scales of our sim- 
ulations. Here, in the majority of the simulation runs, we do 
not observe the formation of a crystal. 




Fig. 5 Two sequences of snapshots of the crystallization process. 
P* = 27 (upper row) and P* = 40 (lower row). Ellipsoids with 
w/, < 8 are shown as dots only. We zoom in to show the relevant part 
of the system. Colors indicate the direction of the principal axis 
(online). 

FIG. [5] shows typical sequences of snapshots. For P* = 27, 
the crystal has a compact structure and nucleation is a rare 
event. For P* = 40, the crystal phase forms a percolating net- 
work. Hence despite the approach to glassy dynamics, we 
observe the classical extremes of kinetics at a first order phase 
transition: nucleation and spinodal decomposition. (A similar 
observation has also been made for hard spheres—.) 
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Fig. 6 (a) Nucleation rate densities for hard spheres and ellipsoids 
as a function of the chemical potential difference |AjU| between the 
meta-stable melt and the stable crystal phase. To allow for direct 
comparison with simulations with other types of local dynamics, the 
rate densities are multiplied by /D, where D is the long-time 
self-diffusion constant. Data is averaged over 25 independent 
simulation runs and errorbars are calculated as the standard 
deviation. 



For P* < 30, where the system crystallizes via nucleation, 
we have computed nucleation rate densities (see FIG. [6]) and 
compared them to values for suspensions of hard spheres (the 
hard sphere data was taken from refs.i^^*^). At equal over- 
compression, the nucleation rate densities for the two systems 
coincide within the errorbars. We conclude that, at these mod- 
erate spect ratios, the orientational degree of freedom of the 
ellipsoids does neither have a strong influence on the interfa- 
cial tension nor on the crystallization process. 



3.3 Dynamic Heterogeneities 

In the rest of this paper, we focus on the influence of dynamic 
heterogeneities on the initial stage of the crystallization pro- 
cess. From each trajectory, we picked the first cluster of ellip- 
soids with nh> 5 that reached a size of 500 particles and stud- 
ied the properties of these particles backwards in time. The 
question we would like to address here is whether we can see 
any difference in the dynamical structure of these particles as 
compared to the rest of the system just before crystallization 
sets in. (We relaxed the criterion for crystallinity here from 
rih > ^ to rih > 5 in order to take particles on the surface of a 
cluster into account. This allows us to analyse our data also 
with resepct to the surface effect discussed in ref. ^'^) 

FIG|7] shows the probability distribution of the distance 
traveled by a particle during a time interval of 5 • 10^ MCsteps, 
i.e. its "mobility". Data for the surrounding liquid is shown as 
circles. The squares present the mobilities of those particles 
that are going to crystallize, taken just before crystallization 
set in, i.e. when 5% of the particles had > 8. There is no 
difference between the two distributions that would allow for 
an identification of the region that is going to crystallize. The 
hypothesis that the most mobile particles crystallize first, be- 
cause they sample their phase space most rapidly, is not sup- 
ported by our data. 

For hard spheres enhanced mobility at the surface of the 
crystalline cluster has been reported and interpreted as a 
signature of a feed back process due to the freeing up of vol- 
ume^. To test for this effect, we computed the mobility when 
80% of the particles had crystallized, but only took particles at 
the surface of the crystallite into account (i.e. with < 8). In 
contrast to refs.^'^, for low pressures we find a shift to lower 
mobilities, and at higher pressures no shift at all. The same 
behaviour is observed if we study the length of the trajectory 
segment instead of the short- time mean- square displacement. 

This observation is confirmed by an analysis of the single 
particle free volume via Voronoi decomposition, see FIG. [8] 
(For a computational routine see ref.— ). We observe no in- 
crease in the single particle free volume at the interface of the 
crystallite. There is no evidence of the modified crystallization 
process described in ref.-. 

Next we look for cooperative motion. We draw a vector 
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Fig. 7 Mobility in the surrounding liquid (circles) in comparison to 
the mobility of particles that are going to crystallize (squares) and of 
the particles at the surface of the crystal, once the crystallite has 
formed (crosses). 



Fig. 9 Alignment of motion of neighbouring particles in the 
surrounding liquid (circles) in comparison to particles that are going 
to crystallize (squares) and particles at the surface of the crystal, 
after the crystallite has formed (crosses). See main text for definition 
of 0. 
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Fig. 8 Single particle free volume in the surrounding liquid (circles) 
in comparison to the single particle free volume of particles that are 
going to crystallize (squares) and of the particles at the surface of 
the crystal, after the crystallite has formed (crosses). 



Ar(/) between the first and the last position of a particle tra- 
jectory travelled over an interval of 5 • 10^ MCsteps. 6 denotes 
the angle between Ar(/) and Ar(j) for neighbouring particles 
/ and j (i.e. if cos = 1 the particles have moved in paral- 
lel.) FIG. [9] shows the distribution of cos(0) for the surround- 
ing liquid (circles), for particles that are going to crystallize 
(squares), and for the particles at the surface of the crystal, af- 
ter the crystallite has formed (crosses). There is no difference 
between these cases. We conclude, in particular, that particles 
that move in parallel are not more likely to crystallize than 
others. 

To conclude we analyse the structural properties of the 
emerging crystallites by evaluating the absolute value of the 
bond-orientation order-parameter |^6(0I ^^^^ particle iM. 
FIG. [TOl shows the distribution of |^6| for particles in the sur- 



rounding liquid (circles), particles that are going to crystallize 
(squares) and of the particles at the surface of the crystal, af- 
ter the crystallite has formed (crosses). Clearly the crystalline 
particles have a higher value of |^6(0I' t>ut there is no differ- 
ence between the particles that are going to crystallize and the 
surrounding liquid. 
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Fig. 10 Distribution of values of particles in the surrounding 
liquid (circles), of particles that are going to crystallize (squares) 
and of the particles at the surface of the crystal, after the crystallite 
has formed (crosses). 



4 Conclusion 

To summarize, we have studied crystallization in suspensions 
of hard ellipsoids at a moderate aspect ratio, a/b= 1.25. We 
picked this system, because it is simple and monodisperse, its 
phase diagram is well known, it shows a glass transition and 
it can be realized experimentally. We pressure-quenched the 
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system into a meta- stable melt state and studied its dynamical 
and structural properties during the subsequent crystallization 
process. For high overcompression, we observe two step re- 
laxation, sub-diffusive behavior, and increasing dynamic het- 
erogeneities. For moderate amplitudes of overcompression 
P* < 30, the system remained in a meta-stable melt state for 
long times and then crystallized via nucleation and growth. 
We showed that the nucleation rate densities for hard ellip- 
soids are the same as for monodisperse hard spheres (within 
the error-bars). For the glassy regime at high overcompres- 
sions, crystallization sets in on a time scale comparable to the 
relaxation times of the dynamic structure factor and particles 
diffuse less than their diameter. The structure of the crystalline 
regions resembles a percolating cluster. Hence, despite the 
glassy dynamics, we observe the classical extreme cases of 
crystallization: nucleation and spinodal decompositon. 

In order to test for conjectured correlations between dy- 
namic heterogneities and crystallization events, we identified 
the regions that crystallized first and followed their behaviour 
backwards in time. We did not find any signature in the dy- 
namic structure of the melt that would allow to predict future 
crystallization sites. 

Ref.^ reported signatures of enhanced mobility at the 
melt/crystal interface for hard sphere glasses due to the freeing 
up of volume. We did not see this effect. For high overcom- 
pression, ellipsoids close to the melt/crystal interface are as 
mobile as ellipsoids in the rest of the melt and for low over- 
compression they even slow down. This observation is con- 
sistent with our analysis of the single particle free volume via 
Voronoi cell decomposition. In addition, we tested for coop- 
erative motion and did not find any signal that would allow to 
predict a crystallization site. 

For our study, we conclude that there are no significant cor- 
relations between glassy dynamics and the crystallization pro- 
cess. 
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